Building integral projection models with nonindependent vital rates

Abstract Population dynamics are functions of several demographic processes including survival, reproduction, somatic growth, and maturation. The rates or probabilities for these processes can vary by time, by location, and by individual. These processes can co‐vary and interact to varying degrees, e.g., an animal can only reproduce when it is in a particular maturation state. Population dynamics models that treat the processes as independent may yield somewhat biased or imprecise parameter estimates, as well as predictions of population abundances or densities. However, commonly used integral projection models (IPMs) typically assume independence across these demographic processes. We examine several approaches for modelling between process dependence in IPMs and include cases where the processes co‐vary as a function of time (temporal variation), co‐vary within each individual (individual heterogeneity), and combinations of these (temporal variation and individual heterogeneity). We compare our methods to conventional IPMs, which treat vital rates independent, using simulations and a case study of Soay sheep (Ovis aries). In particular, our results indicate that correlation between vital rates can moderately affect variability of some population‐level statistics. Therefore, including such dependent structures is generally advisable when fitting IPMs to ascertain whether or not such between vital rate dependencies exist, which in turn can have subsequent impact on population management or life‐history evolution.


| INTRODUC TI ON
Population models use estimated (or assumed) vital rates at the individual level to understand many aspects of a population's ecology and evolution: its long-term abundance trajectory and age, size, or state distribution; its sensitivities and elasticities relevant for management; and its optimal life-history strategy, among others.
Variation in vital rates can have important affects on populations (Hamel et al., 2018;Vindenes & Langangen, 2015). This broad concept encompasses variation across individuals, across cohorts, and/ or through time in ways described more below. In many models, potential variation in multiple vital rates is artificially assumed to be independent.
Correlations between growth, survival, reproduction, and/or other traits can change demographic conclusions .
For example, whereas independent temporal heterogeneity in vital rates has been generally predicted to decrease population growth rate, it can actually increase population growth rate when multiple vital rates are correlated (Doak et al., 2005). A completely different example is that persistent individual heterogeneity in vital rates can reveal different optimal life-history strategies in different environmental conditions (Kentie et al., 2020).
Integral projection models (IPMs) are the framework for discretetime population dynamics with continuous individual-state variables (e.g., mass, size) (Easterling et al., 2000). Compared to age-or stagestructured matrix population models, which track abundance for each discrete state category, IPMs track abundance as a distribution (density) for continuous state values. This enables IPMs to more accurately represent populations in which continuous state variables are important predictors of individual dynamics such as growth, reproduction, and survival (Ellner et al., 2016;Merow et al., 2014;Rees et al., 2014). Thus, it may be important to incorporate both variation in vital rates and correlations among multiple vital rates into IPMs.
To what extent have correlated vital rates been incorporated into both estimation and analysis of IPMs? At a basic level, correlation in individual vital rates arising from stochastic life trajectories is almost inherent to a nontrivial IPM. For example, in a size-structured IPM, correlation in growth and survival will arise when both depend on size and individual size trajectories vary due to stochastic growth.
Temporal correlations among vital rates (e.g., a good year is good for each of growth, survival and reproduction) are captured naturally when year-specific transition kernels are estimated or correlated random effects are estimated (Childs et al., 2004;Hindle et al., 2018;Metcalf et al., 2015). Correlations in individual heterogeneity among multiple traits have been considered for life-history tradeoffs and eco-evolutionary IPMs (Coulson et al., 2021;Kentie et al., 2020).
However, there remains a need for systematic formulation and comparison of multiple kinds of correlated vital rates. This will allow for identification of gaps in statistical estimation and IPM analysis methods and comparison of impacts on demographic conclusions for the same data. Some IPM formulations have been sufficiently general to encompass these kinds of correlations from a mathematical perspective Coulson et al., 2017), but case studies and estimation tools have not been as highly developed.
In this paper, the general concept of nonindependence among vital rates includes three quite different categories: (i) labile individual heterogeneity, (ii) temporal heterogeneity, and (iii) persistent individual heterogeneity. Labile individual heterogeneity refers to differences arising from phenotypic plasticity and the random events of a life course . This is also called dynamic condition (Forsythe et al., 2021) or transient heterogeneity (Brooks et al., 2017). For example, an individual who by luck experiences high-growth conditions in early years may continue to be above average in size throughout its life. Labile heterogeneity can also arise from physiological tradeoffs such as costs of reproduction. For example, if an individual gives birth during the spring, its growth rate over subsequent months may be lower than if it had not given birth. In this example, the heterogeneity could be viewed as an individual-level trade-off between reproducing or growing more, although rigorously proving such causality cannot be done without a controlled experiment (Coulson, 2012;Knops et al., 2007). In statistical models, labile individual heterogeneity can be incorporated by making the transition (projection) kernels for multiple vital rates interdependent. Below, we consider both a standard regression framework and introduce a new copula approach for modelling such interdependence.
Temporal heterogeneity is driven by a shared covariate, which may be observed or unobserved (latent), that affects multiple traits (Compagnoni et al., 2016;Coulson et al., 2011;Hindle et al., 2018;Metcalf et al., 2015;Vindenes et al., 2014). For example, such a covariate could be annual (or breeding-season) food supply that has a positive correlation with both survival probability and fecundity.
Demographic data spanning multiple years would then show a positive correlation between population-level survival and fecundity values. Note that a factor such as food supply could contribute to both temporal heterogeneity-to the extent individuals experience similar growth in a year due to the same conditions-and/or labile heterogeneity-to the extent individuals experience different growth due to heterogenous food conditions in the same year. We will present two different approaches for modelling correlated temporal heterogeneity, one being to explicitly include a shared and measured covariate that affects multiple vital rates and the other being to implicitly include shared but unmeasured covariates by including correlated temporal random effects. Persistent individual heterogeneity in multiple traits refers to between-individual differences that last their entire life (Brooks et al., 2017). This is also called fixed condition (Forsythe et al., 2021) or heterogeneity (Steiner et al., 2010). For example, one individual's average growth and fecundity rates could remain consistently higher than another individual's rates due to fixed heterogeneity. Persistent individual heterogeneity can be as simple as an univariate quality affecting a single trait (Ellner & Rees, 2006) or as complicated as a multivariate vector affecting the duration of the different life stages of an individual (de Valpine et al., 2014). Persistent individual heterogeneity is necessary to represent genetic variation in models of ecoevolutionary dynamics Vindenes & Langangen, 2015), but it can also represent only phenotypic variation potentially shaped by good site conditions at birth, for example. Processes such as energy acquisition allocation (van Noordwijk & de Jong, 1986) or reproductive strategy trade-offs (Benton & Grant, 1999) could be considered as labile heterogeneity and/or persistent heterogeneity in different cases. In this paper, the statistical models of correlated persistent individual heterogeneity use correlated individual random effects (Brooks et al., 2017;Knape et al., 2011), although they can also use individual-level covariates (Moyes et al., 2011). In summary, the three kinds of individual heterogeneity are biologically and statistically distinct, at least in principle.
Numerous IPM studies have incorporated one or more type of heterogeneity in vital rates, but few have incorporated nonindependent forms of heterogeneity (beyond the correlated vital rates arising from a basic IPM formulation). For example, Ellner and Rees (2006) incorporated persistent and labile individual heterogeneity without correlation, and Ellner and Rees (2006) incorporated temporal heterogeneity without correlation. As described by Vindenes and Langangen (2015), some studies include heterogeneity in estimation but then use only mean traits for analysis and prediction.
Evolutionarily explicit IPMs have included both quantitative genetic traits and phenotypes as state variables, which together can be a kind of correlated persistent heterogeneity Coulson et al., 2017Coulson et al., , 2021Rees & Ellner, 2019). Although these have mathematical similarity in IPM formulation, they are distinct in goals and statistical parameterisation methods compared with a nonevolutionary model with correlated individual traits. Kentie et al. (2020) considered correlated persistent heterogeneity among growth, survival, and reproduction, although they did not estimate these in a hierarchical statistical modeling framework as we do here. It is important to realize that each kind of correlated heterogeneity introduces different implementation challenges both for estimation and for IPM analysis involving multidimensional numerical integration, discussed more below.
Statistical estimation of different forms of nonindependent vital rates can draw on methods from other kinds of ecological analyses that, in some cases, have not typically been used for parameterization of IPMs. For labile individual heterogeneity, one current phenotypic value can be used to predict changes in another, which is basic to the formulation of IPMs. Such dependence can in principle include time lags, although these are not explored here. A potential limitation of the simple regression approach is that correlation among vital rates can be induced only be modifying the marginal distribution of the traits. We introduce the use of statistical copulas in this context as an alternative way to model labile correlations. For correlated temporal heterogeneity, one can include correlated temporal random effects or shared explanatory variables (Evans & Holsinger, 2012;Hindle et al., 2018;Metcalf et al., 2015). Alternatively, one can estimate different kernels for each of many years (Childs et al., 2004).
Relevant to persistent individual heterogeneity, statistical models for individual demographic data routinely include random effects for individual heterogeneity, and multivariate random effects can be correlated ( van Bonnet & Postma, 2016;de Pol & Verhulst, 2006). In the case of marked animals with imperfect detection or recapture, capture-mark-recapture methods can also incorporate correlated individual random effects (Cam et al., 2013;Gimenez et al., 2018).
In this paper, we systematically present statistical methods to estimate different kinds of correlations in vital rates and incorporate those correlations into IPMs. We give methods for modelling correlations in vital rate arising in each of the three categories of heterogeneity, including a new copula method for individual heterogeneity. We show how the methods can be used in a hierarchical statistical framework and discuss some of the computational and implementation challenges involved. In a case study with Soay sheep data, we illustrate that the same data can imply different demographic conclusions when different kinds of correlated vital rates are considered. In addition, even when including correlations does not change point results such as population growth rate or elasticities, it can change the width of uncertainty (credible or confidence interval) propagated from uncertainties in parameter estimates.
The structure of this paper is the following. We begin with a general description of IPMs (Section 2.1) and consider IPMs with independent vital rates (Section 2.2). We next discuss the area of primary focus: IPMs with heterogeneous and nonindependent vital rates (Section 2.3). We note here that while dependency and correlation are not exactly equivalent, we will use the terms interchangeably because of common practice. This is followed by a description of simulation studies and a case study using data from a population of Soay sheep (Ovis aries) in Scotland (Sections 2.5 and 2.6). The results of these studies (Section 3) focus on differences arising from the nonindependent vital rate models on (i) the log population growth rate and (ii) population growth rate elasticities. We conclude with a discussion of the implications of the proposed methods (Section 4).

| General integral projection models
We begin with a description of a family of IPMs that permit the incorporation of temporal, persistent and/or labile individual heterogeneity, using the notation from Childs et al. (2016). Let x denotes the individual state variables, hereafter called "i-states." The i-states comprise labile traits that vary over the life cycle in response to the environment such as body mass, length or breeding status (Coulson, 2012;Merow et al., 2014;Rees et al., 2014). In addition, individuals are further characterised by "q-states," denoted by z. The q-states comprise unmeasured, nonlabile characteristics that are fixed during the lifetime of the individual. In this article, we assume that (i) individuals can be uniquely characterized by (x, z), which essentially assumes that individuals with the same (x, z) are interchangeable, (ii) all vital rate models depend on x, and (iii) selected vital rate models depend on z. The values of (x, z) at one discrete time step later are denoted as (x′, z′).
The state of the population is described by the abundance density, denoted n(x, z, t). The abundance density is defined such that the number of individuals at time t with states in a small interval (x, z) to (x + Δx, z + Δz) is approximately n(x, z, t)ΔxΔz. The total abundance at t can then be expressed as N t , such that The projection of the abundance density over time is described by the following equation, is the time-varying projection (transition) kernel, i.e. the density of individuals evolving from (x, z) to (x′, z′) (Ellner & Rees, 2007). The term d t denotes measured and/or unmeasured timespecific environmental conditions that account for temporal variation. are singletons, and for simplicity, we ignore twinning (Coulson, 2012).
In the following sections, we discuss different ways to construct vital rate models when rates are independent or dependent, given the i-states, x. Motivated by reproduction cost (Gittleman & Thompson, 1988;Tavecchia et al., 2005), we restrict attention to the dependence between growth and reproduction.

| Independent vital rate models
Before describing different formulations of vital rate models, we introduce some additional notation. To begin, we assume that there is only one element in the labile traits, x, and that is the natural logarithm of body mass. For individual j at time t, let m j,t denotes the log body mass (given survival); a j,t the alive (1) vs dead (0) state; r j,t the reproductive (1) vs nonreproductive (0) state (given survival); and c j,t the offspring log body mass (given reproduction). The discrete times are t = 1, ⋯, T.
In terms of parameters, fixed-effect parameters are referenced as with subscripts defining the vital rate and the variable they influence, respectively. For instance, g,0 is the intercept for the growth model and s,m is the slope for the survival model corresponding to the variable m. Also, residual (nonrandom effect) variances are denoted by 2 with the subscript defining the vital rate. In addition to fixed effects, we consider random effects on year and individual for temporal and persistent individual heterogeneity, respectively. These random effects are placed on the growth and reproduction models to capture the potential dependence of interest.
The unobserved temporal or individual random effects are denoted by u and v, respectively. For example, u b,t is the reproduction random year effect in year t, while v g,j is the growth random individual effect on individual j. Random effect variances are denoted by 2 and 2 , and correlation parameters by and , respectively.
Assuming independence between vital rates, parameters for each vital rate model can be estimated separately. For that case, we summarize three of the most commonly used approaches to formulate vital rate models.

| Vanilla Model (I1)
We initially define the "vanilla model," denoted as model I1, as the widely used approach where the vital rates depend only on the labile phenotype, x, corresponding to the log body mass (m) in our Soay sheep example (Easterling et al., 2000;Ellner & Rees, 2006). In particular, parameters are estimated given the individual-level demographic data such that, is the inverse of the logistic transformation. To apply the vanilla model to the projection kernel in Equation (3), we rearrange the vital rate models such that, where ϕ(a; , 2 ) denotes the density function of N( , 2 ) evaluated at a. Here, x = m, and there is no z or d t . The equation for h( ⋅ ) represents an inheritance or the "parent-offspring phenotypic similarity" function (Coulson et al., 2021), with offspring size depending on parent size. For the following models, we assume the same vital rate models as described above if they are not mentioned in the model description.

| Temporal Heterogeneity (I2)
Models with temporal heterogeneity connect vital rates with timevarying factors, such as resource availability, natural enemies, and abiotic conditions. We consider a hierarchical model with independent random effects (Bolker et al., 2009;McCulloch & Searle, 2001) such that, where the random effects u b,t and u g,t are independent to avoid inducing dependence between different vital rate models.
Similar to Equation (5), the vital rate models are rearranged such that, , and there is no z.

| Persistent Individual Heterogeneity (I3)
The persistent individual heterogeneity model, denoted I3, differs from the temporal heterogeneity model (I2)  where the random effect distributions are independent to avoid inducing dependence. In this case, the vital rate models are re-arranged as, where v o b and v o g denote the random individual effects for the offspring. Here, x = m, z = (v b , v g ), and there is no d t . We assume offspring size depends on parent size while offspring random effects are independent of parent random effects.

| Nonindependent vital rate models
We now discuss different ways to induce the dependence structure between vital rate models. Corresponding to the three types of heterogeneity are three categories of models, with a category representing labile individual heterogeneity having two models (D1a and D1b), the temporal heterogeneity category having two models (D2a and D2b), and the persistent individual heterogeneity category having one model (D3).

| Labile Individual Heterogeneity (D1a and D1b)
Models in this category extend the vanilla model I1 to create dependence between reproduction and growth. We construct two types of dependent vital rate models: (i) the reproduction conditional model and (ii) the copula model. The former model treats breeding status as a covariate within the growth model, while the latter model utilizes the copula structure to jointly model growth and reproduction. The latter necessitates estimating multiple kernel functions together, while the former does not.

D1a. Reproduction conditional model
This approach models the growth rate of an individual as a function of the breeding status. In particular, the binary variable, r t+1,j , is a covariate in the growth model such that, Integrating out r j,t+1 to obtain the marginal growth model for the projection kernel, we note that, where the marginal growth distribution is now a mixture of two Gaussian distributions and hence potentially bimodal. Here, x = (m, r), and there is no z and d t .
This model induces a dependency between growth and reproduction that is reflected in the covariance, . This covariance is maximized when b(m) = 0.5 and minimized as b(m) approaches 0 or 1.

D1b. Copula model
Copula methods are a popular approach to construct a joint distribution for correlated random variables given assumed marginal distributions [see, e.g., Chapter 6 of Song, 2007). These models extend univariate linear models to general multivariate models with vector responses and provide a flexible approach to the regression analysis of correlated discrete, continuous, or mixed responses (Anderson et al., 2019;de Valpine et al., 2014).
The copula method relies on Sklar's theorem (Sklar, 1959) which states that any multivariate distribution can be constructed by combining the marginal distributions with a suitable copula function describing the association between the variables. Mathematically, given of variables Y 1 , ⋯, Y n , and a copula function C, the joint CDF can be expressed as, There are a variety of copula functions available that permit different behaviours of multidimensional distributions and typically lead to different dependence structures. However, the marginal distributions of the random variables remain the same irrespective of the choice of copula function. We use the Gaussian copula function to handle the dependence structure for simplicity (Nelsen, 2006;Song et al., 2009). The Gaussian copula function is defined such that, m j,t+1 |m j,t , r j,t+1 : N( g,0 + g,m m j,t + g|r r j,t+1 , 2 g ).
F 1,⋯,n (y 1 , ⋯, y n ) = P(Y 1 ≤ y 1 , ⋯, Y n ≤ y n ) = C(P(Y 1 ≤ y 1 ), ⋯, P(Y n ≤ y n )), where Φ −1 ( ⋅ ) denotes the inverse CDF of a standard Gaussian distribution; Φ D ( ⋅ ) and ϕ D ( ⋅ ) are the CDF and density, respectively, of a n-dimensional Gaussian distribution with a zero vector as mean and covariance matrix D. The diagonal elements of D are all scaled to unity without the loss of generality.
As an example, we briefly describe the copula model used in the Soay sheep case study for correlated growth and reproduction, involving the combination of a continuous and discrete random variable. In particular, we use the Gaussian copula function with a normally distributed random variable for growth, Y 1 , and a Bernoullidistributed random variable for reproduction, denoted Y 2 . Note that the density function and CDF of Y 1 are expressed as, where is the expected value of Y 1 , and 2 is the variance of Y 1 . For the reproduction (Bernoulli) variable, as the raw scale is discrete, we introduce an auxiliary variable X, which is distributed as an uniform distribution (i.e., X: U[0, 1]) and define the new random variable Y 3 = Y 2 + X .
The probability mass function for Y 2 , the probability density function for Y 3 , and the CDFs for both are then expressed as, where q = Pr Y 2 = 0 . Combining Equations (13) and (15), we derive the joint density of Y 1 , Y 3 such that, We can then substitute the growth and reproduction model for Y 1 and Y 2 to obtain their corresponding joint density for parameter estimation. The notation becomes x = (m, r), and there is no z and d t .
Despite the appealing features of copula models, IPMs with copula models give the same projection kernel as the vanilla model, which leads to the identical projection of the population dynamics. This is true because (i) correlations in the copula model do not modify the marginal distributions and (ii) the involved vital rate models (reproduction and growth) are an additive structure. Further details are presented in Appendix S1. Demographically, population change is the same whether individuals who grow less are the ones who reproduced more or not.
However, as discussed more below, the copula remains interesting because it may give different answers for life-history questions involving trade-offs, or estimated parameters may be different, or it may give different kernels when used with time lags or other extensions.

| Temporal heterogeneity (D2a and D2b)
These models induce dependence on vital rates by the time-varying factors, extending the independent temporal heteroegeneity model,

I2.
In particular, when the conditions of a given year are "good" for both growth and reproduction, temporal heterogeneity will create positive temporal correlation among these vital rates, which may generally be the case (Hindle et al., 2018). We consider two models: (i) the shared

D2a. Shared drivers model
This approach includes observed time-varying covariates in the regression functions for vital rate models (van Benthem et al., 2017;Dalgleish et al., 2011;Simmonds & Coulson, 2015). Common choices include environmental indices, e.g., North Atlantic Oscillation, precipitation, temperature, etc. To quantify the additional influence of the drivers on the vital rates, let q t denotes the vector of covariates with an associated vector of regression coefficients ⋅,q , namely The vital rate models are re-arranged for the projection kernel such that, Here, x = m, d t = q t , and there is no z.

D2b. Correlated random year effect model
The second model extends the independent temporal random effects model (model I2). Generalizing these hierarchical models by allowing for dependencies in the random effect distributions induces dependencies between vital rates (Hindle et al., 2018;Metcalf et al., 2015) such that, f 1 (y 1 ) = ϕ(y 1 ; , 2 ) (17) r j,t+1 |m j,t , q t : Bernoulli logit −1 b,0 + b,m m j,t + b,q q t m j,t+1 |m j,t , q t : N g,0 + g,m m j,t + g,q q t , 2 g .
The vital rate models are re-arranged for the projection kernel such that, Here, x = m, d t = u b,t , u g,t , and there is no z.

| Persistent individual heterogeneity (D3)
Similar to the temporal heterogeneity, the model in this category extends model I3 to induce dependence between vital rates for the persistent individual heterogeneity case.

D3. Correlated random individual effect model
We consider a hierarchical model with dependent random effects distribution, similar to model D2b. In particular, we specify, The vital rate models are re-arranged for the projection kernel such that, where ind ( ⋅ ) is the density function of the random individual effects distribution and specified in the last part of Equation (21). Here, x = m, z = v b , v g , and there is no d t .

| Comparison of the models
In Figure 1, we present a graphical representation of the differences between the proposed heterogeneity models. In each of the four scenarios, the individual growth model, g ( ⋅ ), depends on exactly one factor.

| Hybrid models
The proposed models can occur individually or be combined within and/or between the categories (labile individual, temporal, and Prior distributions for all parameters are set to be noninformative and are presented in Appendix S2. We use the trace plot and Brooks-Gelman-Rubin statistic to assess convergence (Gelman & Shirley, 2011). Chains with a value of Brooks-Gelman-Rubin statistic being less than 1.05 are treated as converged.

| Approximation of log s
We use the asymptotic log population growth rate, log , as one metric to compare models. Mathematically, is defined as where N t is the population abundance and can be approximated by solving the integral in Equation (2). It has been shown that log converges asymptotically, even in the temporally stochastic case (Ellner & Rees, 2007).
The log population growth rate of IPMs without temporal heterogeneity can be approximated via the midpoint rule (Easterling et al., 2000). To briefly illustrate the mid-point rule, the projection kernel is discretized into a projection matrix by a sufficient number of mesh points that are of uniform length to discretize (x, z) (Ellner & Rees, 2006). The population growth rate is then obtained as the leading eigenvalue of the projection matrix (Caswell, 2001). Alternatively, we can consider using mesh points that are uniform quantiles of z as the distribution of z is known.
However, when the IPMs include temporal heterogeneity, the midpoint rule becomes inapplicable. In this case, we use the simulation technique of "element-selection" to approximate the log population growth rate (Ellner & Rees, 2007;Rees & Ellner, 2009).
This approach creates a series of projection matrices, K t with the population abundance N t obtained by repeatedly multiplying the projection matrices with a discrete approximation of n (x, z, t). The (stochastic) log population growth rate is approximated using the empirical mean given by, where data in the first L 0 < L years are excluded as transient dynamic to reduce the influence of random initialization. We note that this estimator carries an extra variability caused by finite simulation. Ellner and Rees (2007) showed that the estimator converges to a normal distribution such that, In addition to the log s itself, we are also interested in the variability on log s caused by parameter uncertainty. This parameter uncertainty can be easily propagated within the Bayesian framework since we are able to obtain samples from the posterior distribution of the parameters, which in turn can be used to calculate the value of log and hence obtain summary statistics of the posterior distribution.

| Sensitivity and elasticity analysis
We also estimate the sensitivity and elasticity of the asymptotic log growth rate, log s , with respect to selected vital rate parameters (Rees & Ellner, 2009;Tuljapurkar, 1990;Vindenes et al., 2014). In particular, we note that Coulson et al. (2005) suggests that models incorporating between-process correlations may alter the sensitivity estimate, which in turn has implication for management decisions.
Here, we apply a central-differencing approach to approximate the sensitivity such that, where s ( + ) is the estimate of s when the target parameter equals to + . By running preliminary tests, we found that = 0.005 is small enough to give precise estimate for all sensitivities of interest. Given the estimate of sensitivity, elasticity of is obtained as, We note that the sensitivities/elasticities of the copula model (D1b) are the same as for the vanilla model (I1), similar to . To see of Ellner et al., 2016) such that, where both terms in the integral remain unchanged because the copula model does not distort the marginal vital rate models.

| Simulation study
We conducted a simulation study to investigate how sensitive the  (Table 1) or fixed (Table 2).
Randomly generated parameters allowed us to look at how summary statistics change over small ranges of variation in a coarse way, without looking at changes in relation to each parameter one by one.
The distributions and values are motivated from the data in the case study, but slightly adjusted to show the difference between models with and without correlations.
The simulation study looks at theoretical behavior of the IPM models, not at statistical properties of parameter estimation. It reveals how model summary statistics shift with particular parameters, but not how parameter estimation performs if the wrong model is fitted to the data. Within the simulation study, we compare the independent models (I1 − I3) and three of the dependent models (D2a, D2b, D3). We do not include the models with labile individual heterogeneity as: (i) the impacts on log by the reproduction conditional models (D1a) are always negative when ′ < 0, and (ii) the copula model (D1b) and vanilla model (I1) are theoretically equivalent due to the unchanged marginal property (given the same parameter values). For models with temporal heterogeneity, we set L 0 = 1000 and L = 10, 000.

| Soay sheep case study
We apply the different models to data on Soay sheep. The individuallevel demographic data consist of information from marked female sheep in the Village Bay area on the island of Hirta in the St. Kilda archipelago, Scotland, from 1986 to 1996. Details of the Soay sheep and data collection protocol can be found in Clutton-Brock and Pemberton (2004), and the data are available from Coulson (2012).
Using preliminary runs for the estimation of parameters of the vital rate models, we set the burn-in and total iteration numbers for the MCMC algorithm to be 20, 000 and 100, 000 for the ma- proposed models that intend to capture the correlation between reproduction and growth. In this article, we analytically marginalise out the missing data to estimate parameters of interest.

| Simulation study
In Figure 2, we present the pairwise results of the vanilla model ability is more substantial for models with temporal heterogeneity,  (Figure 2a-c). This is in line with the result that although uncorrelated temporal heterogeneity is generally predicted to decrease log s , correlated temporal heterogeneity can increase log s (Doak et al., 2005;Fieberg & Ellner, 2001). Also, the temporal heterogeneity models and persistent individual heterogeneity model cause different impacts on log s . For example, temporal heterogeneity appears to lead to reduced log s ; similarly, increasing the correlation in temporal heterogeneity models leads to a decrease in log s (Figure 2a,b). However, persistent individual heterogeneity models have the reverse effects ( Figure 2c). Finally, we note that the trend on log s against correlation does not translate into that of elasticities. The decreasing trend of the temporal heterogeneity disappears (Figure 2a

| Case study on Soay sheep
In Appendix S3, we present the posterior summary estimates of the model parameters for different models. Three dependent models

| Comparison of log s
We use 500 parameter values sampled from the posterior distribution to approximate the (stochastic) log population growth rate. The uncertainty from parameter estimation are hence propagated into the posterior distribution of log s . In the temporally stochastic models, we set L 0 = 1, 000 and L = 10, 000 to approximate log s . Table 3 provides the corresponding summary statistics of log s for each model.
We first observe that the mean of log s ranges approximately from 0.03 to 0.04, which translates into a 3 to 4% annual population growth rate. There is considerably more variability, however, in the uncertainty about log s . In particular, the width of the credible intervals of log s by models with random year effects (I2, D2b) are around 35\% larger than that of the rest of the models. Second, we observe that the uncertainty on log s caused by parameter uncertainty is larger than the bias caused by ignoring the correlation structure.
This is similar to the empirical result of Compagnoni et al. (2016) that parameter uncertainty outweighs the bias caused by ignoring the correlation structure. Further, we note that log of the vanilla model (I1) and the copula models (D1b) are slightly different despite the theoretical equivalence between the IPMs. This is because the parameter estimates between the models are different.
Finally, we note that the predictions of the shared drivers IPM (D2a) depend on the distribution of the winter NAO. Adjusting the distribution of the winter NAO may lead to different distributions of log s hence interpretation. In Appendix S4, we consider three other distributions obtained by using a nonparametric bootstrapping approach of the NAO in different years.

| Comparison of elasticity
We approximate the elasticities of four parameters, again using the sampled parameter values from the posterior distribution, presented in Table 4. We observe that models with random temporal effects lead to a larger variability in the elasticities, which is similar to the observation in log s . Additionally, we note that the correlated random in- Unlike independent IPMs, these modelling approaches explicitly characterise the dependency between vital rates, permitting the quantification of between-process correlation. As a data-driven method, this is better than assuming either no correlation, or perfect correlation across vital rates, i.e., assuming the correlation coefficient to be 1 or − 1 (Benton & Grant, 1999;Coulson et al., 2011).
Amongst the proposed methods, application of the copula method for modelling vital rates is novel to IPMs. However, given the same estimates for the common parameters, the dependence structure of an IPM using copula models may lead to theoretically equivalent projections as the independent (vanilla) IPM. This is because (i) correlations in the copula model do not modify the marginal distributions and (ii) the involved vital rate models (reproduction and growth in our analysis) have an additive structure. In practice, however, copula IPMs will still differ from the vanilla IPMs due to differences in parameter estimates. Further, such theoretical equivalence will not remain with alternative copula structures, for example, when we consider the previous breeding status r j,t as opposed to the current breeding status r j,t+1 in the copula structure with the growth vital rate. It may be appropriate to condition on reproduction at time t for some species, particularly when multiple reproduction-related activities can cause energy loss in the parents including mating, gestation, parturition, lactation, etc (Gittleman & Thompson, 1988).
Also, copula models can be applied to other aspects of IPMs. For instance, the multidimensional random effect distribution can be constructed by copula models, which bring extra flexibility to the models. The use of copula models within this general context is an area of current research.

| Simulation and case study
In the case study of Soay sheep, the different IPM structures yielded relatively similar population estimates. This is most likely because the parameter uncertainty (which was ignored in the simulation studies) outweighed the impact of between-process correlation (Compagnoni et al., 2016). In contrast, the results for both the simulation and the case study show that (i) different models for dependence between vital rates can yield similar (nearly identical) log s , but different elasticities and (ii) variability of the population statistics are moderately affected by the correlation between vital rates.
Random effect models are commonly used to model dependence structures (Dingemanse & Dochtermann, 2013;Vindenes et al., 2014). Based on the simulation study, it appears that temporal and persistent heterogeneity can lead to differences in the estimated target statistics and their associated variability. Results suggest that the variability increases as the correlation increases. This aligns with the general understanding that extreme values are more likely to be generated, and hence, the variability of the target statistics increases when the correlation is large and positive (Doak et al., 2005;Fieberg & Ellner, 2001). Empirical results about the correlation in temporal variation have been discussed previously (Hindle et al., 2018;Metcalf et al., 2015). Additional random effects models can also be investigated, given available data, for example, allowing for nested spatial heterogeneity (Olsen et al., 2016), or independent/ crossed structure of spatial and temporal heterogeneity (Jacquemyn et al., 2010). Such heterogeneity structures can provide additional flexibility and more complicated correlations in vital rates and hence IPMs.

| Recommendation
In practice, model selection procedures are often carried out to determine whether one model is preferable to all others. However, we note that some of the proposed methods (D1a, D1b) do not allow unbalanced data, whereas other proposed methods (D2a, D2b, D3) are flexible for unbalanced/balanced data (Verbeke et al., 2014). Such differences complicate model selection, which usually assumes the competing models use the exactly same data. This is an area for future research.
In general, incorporating these five (biologically/statistically) distinct methods (in hybrid/separately) in IPMs may be beneficial.
Although the correlations have little impacts on some statistics of interest (e.g., log s ), our empirical results show that elasticities of the unknown parameters and the associated variability are moderately affected by these correlations. These results may provide insights into the relationship between the possible dependencies on individual-level vital rates and target population statistics.
Therefore, we conclude that including such dependent structures is generally advisable when fitting IPMs to ascertain whether or not such between vital rate dependencies exist, which in turn can have subsequent impact on population management or life-history evolution.

ACK N OWLED G EM ENTS
We are grateful to Adam Butler for the helpful discussion. YLF was funded by Biomathematics and Statistics Scotland and University of Edinburgh PhD studentship. RK was supported by the Leverhulme research fellowship RF-2019-299.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest. and writing-review and editing (lead).

DATA AVA I L A B I L I T Y S TAT E M E N T
The demographic data that support the findings of this study are openly available at http://www.oikos journ al.org/sites/ oikos journ al.org/files/ appen dix/sheep_data_1986_to_1996.csv from Coulson, 2012. The NAO data that support the findings of this study are openly available at https://cruda ta.uea.ac.uk/cru/data/nao/nao.dat.
The example code that support the findings of this study are openly available at https://github.com/Eddie Fung/IPM-non-indep enden tvital -rates.